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We study the effects of substructure in the Galactic halo on direct detection of dark matter, on 
searches for energetic neutrinos from WIMP annihilation in the Sun and Earth, and on the enhance- 
ment in the WIMP annihilation rate in the halo. Our central result is a probability distribution 
function (PDF) P(p) for the local dark-matter density. This distribution must be taken into account 
when using null dark-matter searches to constrain the properties of dark-matter candidates. We 
take two approaches to calculating the PDF. The first is an analytic model that capitalizes on the 
scale-invariant nature of the structure-formation hierarchy in order to address early stages in the 
hierarchy (very small scales; high densities). Our second approach uses simulation-inspired results 
to describe the PDF that arises from lower-density larger-scale substructures which formed in more 
recent stages in the merger hierarchy. The distributions are skew positive, and they peak at densities 
lower than the mean density. The local dark-matter density may be as small as l/10th the canonical 
value of ~ 0.4 GeV cm -3 , but it is probably no less than 0.2 GeV cm -3 . 

PACS numbers: 95.35.+d,98.35.Gi, 98.35. Pr, 98.62.Gq 
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I. INTRODUCTION 

The nature of dark matter remains a mystery. Among 
the plethora of dark-matter candidates, there are two 
classes that are sufficiently promising to motivate major 
experimental searches. The first is a weakly interact- 
ing massive particles (WIMP) 0,0, which may arise in 
supersymmctric extensions of the standard model or in 
theories that include universal extra dimensions (UED) 
3 . The other leading candidate is the axion [1] , hypoth- 
esized to solve the strong-CP problem. 

WIMPs may be detected indirectly through observa- 
tion of gamma rays or cosmic-ray positrons, antiprotons, 
and/or antideuterons produced when WIMPs annihilate 
in the Galactic halo or in the halos of some extragalac- 
tic systems. WIMPs might also be detected via obser- 
vation of energetic neutrinos produced by annihilation 
of WIMPs that have accumulated in the Sun and/or 
Earth. There is also an effort to detect dark-matter par- 
ticles in low-background experiments via detection of the 
0(100 keV) energy they impart to nuclei from which they 
clastically scatter. Likewise, dark-matter axions are be- 
ing sought directly by conversion to photons in resonant- 
cavity experiments [5J. If dark matter is composed of 
WIMPs, we may also get some clue to its nature from 
forthcoming LHC accelerator experiments @ . 

Predictions for the event rates for any of these detec- 
tion schemes depend on the distribution of dark matter 
in the Galactic halo. The flux of gamma and cosmic 
rays from WIMP annihilation depends on an integral of 
the square of the dark matter density over volume, while 
the rates for energetic neutrinos and direct detection de- 
pend only on the local dark-matter density. In early cal- 
culations, it was assumed simply that dark matter was 
smoothly distributed in the Galactic halo, with either an 
isothermal or NFW radial density profile, and with a 



local dark-matter density fixed by Milky Way dynamics 
to be pq ~ 0.4 GeV cm -3 . 

It has become clear in recent years, from theory and 
N-body simulations, that the distribution of dark matter 
in the Galactic halo is not likely to be perfectly smooth, 
and that some of the dark matter in the Galactic halo 
will be distributed into subhalos with a variety of sizes 
Q . This substructure hierarchy may extend to very small 
scales H, [13 EH El, El] - For WIMPs in supersymmetric 
and UED models, the cutoff scale of the power spectrum 
is in the range [10 -6 - 10 2 ]M e pj, For axions, the 
growth of perturbations is suppressed on scales smaller 
than (m Q -ff) -1 / 2 [H, [13], the geometric mean of the in- 
verse of the axion mass and the Hubble constant; for 
example, for an axion of mass m a ~ 10 -5 eV, the cutoff 
scale corresponds to about 10 -11 Mq. 

The implications of such substructure for searches for 
cosmic rays from WIMP annihilation in the halo have 
been explored extensively in the literature [l8l |. Simply 
stated, if a dark matter halo contains substructure, then 
the volume integral of the density squared is increased 
by some boost factor B. For example, Ref. [lj| recently 
claimed B ~ 10 3 enhancements in the direction of the 
Galactic anticenter, although Ref. [2(| finds an enhance- 
ment B ~ 2 — 5. An analytic approach presented in 
Ref. [2l| also finds smaller values of B, of order B ~ 10. 

The implications of substructure for direct detection 
of dark matter have, however, been comparatively ne- 
glected. This is a serious omission, as substructure im- 
plies fluctuations in the local dark-matter density, and 
these fluctuations imply an uncertainty in the predicted 
rates for direct-detection experiments and for energetic- 
neutrino searches. A simple but illustrative example 
places all of the dark matter in subhalos of density B 
times the mean density. In this case, the total anni- 
hilation rate is enhanced by a boost factor B, but the 
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probability that the Solar System lives in such a sub- 
halo is only B~ l . If B ;§> 1, this probability is small, 
implying dire prospects for direct-detection, an alarming 
conclusion that warrants further investigation. 

Few past studies focused on the implications of sub- 
structure for the local dark-matter density. Some fo- 
cused on implications derived from numerical simulations 
[22l[23j. thus limited by resolution effects. Others looked 
at the contribution of tidal streams to the local density 
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The goal of this paper is to begin addressing questions 
such as: What are the possible values of the local den- 
sity? How small can they be? What are the most prob- 
able values? To answer these questions, we calculate a 
probability distribution function (PDF) P{p) for the lo- 
cal dark-matter density. This PDF accounts for fluctua- 
tions in the local density due to substructure. The mean 
of this distribution is ~ 0.4 GeV cm~ 3 , the canonical 
value determined from dynamics, but we find that the 
local dark-matter density may be as small as l/10th, but 
probably no less than half this canonical value. The PDF 
we calculate will have to be convolved with the results of 
null searches to infer constraints to particle-dark-matter 
parameters. In addition, the PDF can be used to assess 
what a substructure enhancement in the halo WIMP an- 
nihilation rate implies for the local dark matter density, 
and vice-versa. 

Our first approach, in Section II, uses an ansatz about 
the survival fraction for subhalos that form early at high 
densities. General arguments will be used to bracket the 
plausible range of PDFs, and also the plausible lower 
limit to the local dark-matter density. This analytic ap- 
proach capitalizes on the scale-invariant nature of the 
hierarchical formation of galaxies to extrapolate simu- 
lation results to the highest-density, smallest-scale sub- 
structures that form earliest and that are beyond the 
reach of simulations. We then perform a second calcula- 
tion that assumes substructures with NFW profiles are 
distributed in the Milky Way halo. This second approach 
describes the local PDF due to larger-scale substructures 
which have formed from more recent stages in the merger 
hierarchy, and it thus complements the first approach. 
In Section III, we improve upon this latter calculation 
by using ingredients on subhalo mass functions and con- 
centration parameters taken from numerical simulations. 
We compare our work to past literature in Section IV, 
and we state our conclusions in Section V. 



II. ANALYTIC MODEL 

In hierarchical structure formation, small halos col- 
lapse first, and they then merge to form more massive 
structures. Because each halo virializes to roughly 200 
times the mean cosmological density at collapse, the halo 
density generally increases with decreasing subhalo mass. 

However, a dark matter halo is unlikely to remain com- 
pletely intact as it merges into larger structures in the 



hierarchy, and a significant fraction of the mass of each 
subhalo will be stripped as it is embedded into larger ha- 
los. As a result, a significant fraction of the mass of any 
halo in the hierarchy will be smoothly distributed, rather 
than contained in substructure. Nevertheless, there will 
be some fraction that is contained in subhalos of higher 
density, originating from an earlier stage in the hierarchy, 
and some fraction of those will be entrapped into even 
smaller and denser sub-subhalos, etc., all the way down 
to the smallest scale in the hierarchy. It is important 
to note that N-body simulations are unlikely to ever be 
able to resolve the complete hierarchy of substructures, 
from 10 12 M Q to 10 -10 M@, or even smaller, and so some 
analytic approach is required. Likewise, extended-Press- 
Schechter calculations of substructure keep track of the 
most massive subhalo that each mass element occupies, 
not the least massive subhalo in the hierarchy, which is 
the step in the hierarchy that determines the local den- 
sity. 



A. Local density probability distribution function 

We describe structure formation as a series of steps in 
a hierarchy, where the first stages consist of the lowest- 
mass and densest subhalos, which then merge in subse- 
quent steps in the hierarchy to form more massive struc- 
tures of lower density. Note that this model simplifies by 
assuming that each halo has a uniform density. However, 
as we show below, the qualitative results we obtain are 
essentially unaltered, at least at high density, if we as- 
sume more realistic density profiles, such as isothermal 
or NFW profiles. 

Let /(pi) be the fraction of mass, per logarithmic den- 
sity interval centered on pi, that is not yet locked up in 
halos of higher density. Let F(pi) then be the fraction 
of the mass today in the Milky Way halo that exists at 
density greater than p\. Then F(p\) is related to f(pi) 
through 



dF f , 

— = -J-(l-F). 

dpi pi 



(1) 



In other words, — (dF/dpi)dpi is the fraction of the mass 
in the halo today that has a density in the interval p\ — > 
px + dpi, and this is f / pi times 1 — F, the fraction of 
mass not yet locked up in higher-density subhalos. 

Strictly speaking though, p\ is the density that a 
halo would have if the mass was uniformly distributed 
throughout the halo. The true density p of the smoothly 
distributed component of any subhalo is reduced because 
some of the mass is in subhalos of higher density. A 
given halo has a total mass M and occupies a volume 
V = M/pi, but the mass in subhalos is MF(pi), and 
these subhalos occupy a volume 
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dp'^-dF/dp'JiM/pip'J]. 



(2) 
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The smoothly distributed matter has mass M[l — F(pi)}, 
and it occupies a volume V — V^ u b, and so the true density 
of the smoothly distributed component of a halo of mean 
density p\ is 



p{Pi) = Pi 



I-F(Pi) 



1+ «w sg-dd-X 



(3) 



Thus, for example, locally in the Milky Way, the den- 
sity pq ~ 0.4 GeV cm -3 is the value determined from 
dynamics, but /?, the smoothly-distributed component, 
may be smaller, reduced because some of the mass is in 
subhalos. 

Eq. ([3]) is an integral equation for p(pi). To solve it, 
we differentiate the expression for p(pi) with respect to 
pi to obtain a differential equation, 



dp = P 2 
d Pl pl[l-F( Pl )Y 

This can then be integrated to obtain 



Pi.Pi) 



dp 1 , 



in 



(4) 



(5) 



This result could have also been obtained simply by dis- 
cretizing and summing the contributions of each density 
interval to the volume. 

To recapitulate, we postulate a survival fraction f(pi) 
for halos of mean density p%. We obtain the fraction 
of the Milky Way mass smoothly distributed in halos of 
mean density p\ by solving Eq. ([1]). We then obtain the 
true density p of the smoothly distributed component of 
halos of density p\ from Eq. ([5]). The fraction of mass 
in the Milky Way with density in the interval p — > dp 
is then (-dF/dp) = (-dF/dpi)(dp/dp 1 )- 1 . The final 
step is then to simply note that the density probability 
distribution function P(p) that we seek is actually the 
fraction of the Milky Way volume, rather than mass, at 
density p, and this is 



P(P) 



1 dV 
V~dp 



(l/p)(-dF/dp) 



j^dp'AVpip'M-dF/dp^ 



(6) 



With this result, we have mapped the problem of deter- 
mining P(p) to the problem of determining the survival 
fraction f{pi). In principle, f(p) can be determined with 
simulations, but the appropriate simulations have not yet 
been performed. Thus, here wc will make some educated 
guesses, based partially on estimates from existing simu- 
lations, and then explore the implications of these guesses 
for P(p). 

N-body simulations of cosmological formation of 
Milky- Way-like halos suggest that only about [10-20]% of 
the Milky- Way-halo mass is in its 10 largest progenitors 
[27| . This fraction could potentially be larger as halos 



in future numerical simulations get more highly resolved. 
For example, the survival fraction of early forming sub- 
halos (thus tightly bound) is higher (~ 40% of their mass 
survives [Hj]) than recently forming subhalos; therefore 
if subhalos are resolved at higher densities, the fraction 
of mass in substructure as measured in simulations could 
increase. If we assume [10-40]% of the Milky Way halo 
mass is in subhalos with a density ~ l/10th that of the 
mean density of the Milky Way, then this implies very 
roughly that f{pi) — 0.04 — 0.16, at least for the most 
recent stage in the hierarchy (pi ~ p & ). Therefore, for 
the remainder of this calculation we will consider a range 
/(po) = 0.05 - 0.2. 

How then does f(pi) vary for higher masses? Be- 
fore answering, we first note that the density and mass 
of the inner 10 kpc of the Milky Way, where we live, 
suggest a formation redshift z ~ 10 for the inner 10 
kpc of the Milky Way halo. At these and higher red- 
shifts, the Universe is Einstein-de Sitter. Moreover, 
the effective spectral index for primordial perturbations 
spans the range —2.6 < n e g < —2.1 for mass scales 
10 11 Mq > M > 10~ 6 M ffi , corresponding to the range 
of substructure mass scales we are considering here. If 
we approximate n e g ~ constant, then in the appropri- 
ate range of stages in the hierarchy, structure forma- 
tion is scale invariant. One reasonable guess is thus that 
f(pi) =constant, that the hierarchy of substructures in 
the Galactic halo is scale invariant. 

However, what is more likely, for several reasons, is 
that f(pi) decreases with increasing p\. First of all, 
subhalos in the earlier stages of the hierarchy undergo 
more orbits in their parent halos between the time they 
form and today, and this increases the efficiency with 
which they will be stripped [20|. Secondly, smaller sub- 
halos may be disrupted by tidal interactions with stars 
(28|. Third, the effective spectral index n e s — > —3 for 
earlier stages in the hierarchy. Thus, earlier in the hi- 
erarchy there is less time for a subhalo to virialize fully 
before it becomes enveloped by the next larger halo in 
the hierarchy. To model these effects, we thus consider 
f(pi) = f(PQ)(pi/Pe)- a with a > 0. 

The differential equation in Eq. (fTJ) is readily solved 
for these f(pi) parameterizations. For f(p\) =constant, 
F(p) = 1 — {pi/pmaxY ■ For example, if we take 
p max /p e = 1000, then for f{p x ) = 0.05, 0.1, or 0.2, we 
find that in the Milky Way, F(p e ) = 0.29, 0.5, and 0.75, 
respectively; the remaining (1 — F) mass is in the smooth 
component. Thus, if / = 0.1, then half of the mass of 
the Milky Way halo at radii r ~ 8 kpc is smoothly dis- 
tributed, and the other half is in subhalos. A larger / re- 
duces the smooth component, as does a larger p max . As 
we show below, a power-law f(pi) increases the smooth 
component. The density of the smooth component can 
be estimated to be (1 — F)pQ. More precisely, solution 
of Eq. (UJ) gives the density of the smooth component as 
p(p&) = 0.74, 0.55, and 0.3, respectively. 
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FIG. 1: The local dark-matter-density probability distri- 
bution function P(p) for the analytic model, as a func- 
tion of the density p scaled by the mean density p , for 
{/(pa), «} = {0.05,0} (black solid curve), {0.2,0}, (dotted 
red curve), {0.05, 1}, (short-dash blue curve) {0.2, 1} (long- 
dash green curve). We smooth the Dirac delta function for 
the smooth component to a Gaussian of rms one-tenth the 
smooth-component density. The power-law tails are due to 
subhalos. 



For a power-law /(pi) = f(p®){pi/pe) a (and a > 0), 

/ " 



F{ P i 



1 — cxp < 

a 



1 



Pi 



(7) 



and if a is not too small (a > 0.5 ), the dependence 
on the cutoff density p m ax disappears. In that case, 
F(Pq) ^ 1 - exp[-/(p )/a] ~ f(pe)/a, where the 
last approximation valid for /(p©) <C a. One guess 
for /(pi) is that it is inversely proportional to the for- 
mation time t oc p -1 / 2 for halos of density p. If so, 
then F(po) ~ 2/(p Q ) is the fraction of the Milky Way 
mass that is in subhalos, roughly 20% for /(p©) = 0.1. 
A stronger dependence on pi would result in a smaller 
F(pq). It is then straightforward to calculate the PDF 
P(p) to find power-law dependence P(p) oc p~( 2+Q ) for 
densities above the smooth density. 



B. Results and Discussion 

To illustrate, Fig.rjjshows the PDFs, on a log-log plot, 
for four combinations of the parameters /(p©) and a in 
the model. The distributions feature a high-density tail 
(close to a power-law), due to the fraction (1 — F) of the 
mass that is in substructure. Also shown is a smooth 



component at a density p sm < p©; for purposes of illus- 
tration, we smooth the Dirac delta-function dependence 
of this smooth component to a Gaussian of rms one-tenth 
the smooth-component density. The Figure shows that 
as /(p©) is increased, the amplitude of the high-density 
tail is increased at the expense of the smooth compo- 
nent. We also see that increasing a increases the smooth- 
component density while decreasing the amplitude of the 
high-density tail. 

Table U provides numerical results for the smooth- 
component density, the fraction of the volume occupied 
by the smooth component, and the annihilation enhance- 
ment B (discussed in more detail below) for PDFs P{p) 
parameterized by the survival fraction /(p©) and the 
survival-fraction power-law index a. 

Here are some comments and general conclusions from 
these results so far: 

(1) If /(p©) is roughly / < 0.2 as suggested by nu- 
merical simulations, then the reduction in the smooth- 
component density is no less than 30% the mean den- 
sity. The implied fractional uncertainty in the local dark- 
matter density is thus not too much larger than that 
(roughly factor of two) implied by uncertainties in the 
stellar/gas contribution to the local rotation curve, or 
the uncertainties that arise if we allow for a flattened halo 
(which generally increase the local dark-matter density). 
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TABLE I: The density p S mooth of the smooth component, the 
fraction of the volume occupied by the smooth component, 
and the annihilation enhancement B for density distributions 
P(p) parameterized by the survival-fraction amplitude /(p©) 
and the survival-fraction power-law index a. 



(2) The fraction of the volume occupied by the smooth 
component is larger than the fraction of the mass in the 
smooth component, as the higher-density components oc- 
cupy correspondingly less volume. As a result, in most of 
the models, the density in the vast majority of the halo 
volume is the density of the smooth component. 

(3) We have assumed that each halo and subhalo has 
a uniform density. More realistically, the density of each 
halo and subhalo will decrease with radius, perhaps with 
an NFW profile, which has a density that depends on 
radius r as p oc in the inner regions. The volume in 
such a halo changes with density as (dV/dp) oc p~ 3 for 
p — > co. This falls with density more rapidly at large p 
than P(p) oc p 2+a as long as a < 1. Thus, if a < 1, 
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the high-density scaling of our P(p) will be unaltered 
by convolving our halo density distribution with NFW 
profiles. 

By contrast, this analytic model is weakest perhaps at 
low densities (larger-scale substructures), where details 
of the low-density outer regions of merging substructures 
may be important in determining the local density, the 
density of what we have called the smooth component. 
To address this shortfall, we continue in the next subsec- 
tion to calculate the halo PDF in a different model, and 
one that will then preview the numerical results that we 
present later. 



C. Discrete subhalo populations 

In this Section we consider another illustrative model, 
one in which earlier generations of subhalos are subsumed 
and then smoothly distributed in radius into their larger 
host halos. In this case, the Milky Way halo is seen as an 
amalgamation of smaller halos, distributed in mass with 
some mass function dn/dM (e.g., a Press-Schechter-like 
mass function, or one taken from simulations). Further- 
more, each of these subhalos has an NFW density profile. 
In this scenario, there is no smooth component of dark 
matter distributed in the Milky Way halo; instead, the 
whole halo is made up of adjacent individual subhalos. 
Although there may be a wide variety of halo masses in 
this scenario, they all have comparable mean densities, 
and mean densities comparable to the Milky Way den- 
sity. If this is the case, then the subhalos fill the entire 
volume of the Milky Way. Such a scenario neglects the 
halos-in-halos problem, but it may more accurately de- 
scribe the PDF due to larger-scale substructures, which 
arise from the most recent stages in the merger hierarchy. 

Consider now the volume of radius R occupied by one 
of these subhalos. If the Milky Way mass were uniformly 
distributed, the matter in this volume would have some 
uniform density p. However, the subhalo has its own 
NFW density profile, 



10 r 



p{r) 



(r/r s ){l + r/r a f 



(8) 



where r s is the scale radius. This scale radius is related 
to the concentration parameter c v through R = c v r s ; i.e., 
halos with larger c v are more centrally concentrated. The 
characteristic density p s is related to p, R, and c v through 
p = 3p s f{c v )/cl, where /(cv) = hic v - c v /(l + c v ). 

The fraction of the volume at density p in this model is 
P(p) = (l/V)(dV/dp) = (3/R 3 )r 2 (dr/dp). It can be cal- 
culated analytically, but the expressions are algebraically 
unwieldy and unilluminating. The important thing is 
that (dV/dp) oc p~ 3 for p ^> p s and (dV/dp) oc p~ 2 for 
p <C p s - The results for P(p) in this model are shown 
in Fig. O These models predict a broad range of p, but 
the median densities are p mo d = 0.35 p© and 0.58 p© for 
c v = 10 and 2, respectively. The 95% C.L. lower limits 
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FIG. 2: The local dark-matter-density probability distribu- 
tion function P(p), for the discrete-subhalo model, as a func- 
tion of the density p scaled by the mean density p Q , for 
c v = 10 (black solid curve) and c v — 2 (dotted red curve). 



to p are P95 = 0.20 pq and 0.36 p© for c v = 10 and 2, 
respectively. 

Note that nowhere in this discussion did we specify 
the subhalo mass or mass function. The results apply as 
long as all of the subhalos have the same concentration 
parameter and virial density, independent of their mass. 
More realistically, there will be a range of concentration 
parameters and virial densities, and this will be explored 
in the next Section. Note also that we have assumed 
no smooth component, but more realistically the outer 
parts of each subhalo will be tidally stripped to provide a 
smooth component. This can be seen roughly by noting 
that the density of the NFW profiles at the maximum 
radius is roughly 0.19-0.34 the mean density. This will 
also be explored in the next Section. 



D. Annihilation enhancement due to substructure 

As we mention in the Introduction, substructure in the 
halo implies an enhancement to the dark-matter annihi- 
lation rate. For example, a 10 _4 M© halo has a char- 
acteristic density which is roughly 1000 times the local 
dark-matter density. If all the dark matter in the halo 
today resided in such subhalos, then the annihilation flux 
would be boosted by a factor of 1000. However, the prob- 
ability for the Solar System to be in such a subhalo would 
then be only 0.1%. 

The analytic models developed here allows us to eval- 
uate the enhancement of the annihilation rate due to 
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substructure, and relate it to the local dark-matter den- 
sity. The enhancement in the annihilation rate, over the 
smooth dark-matter distribution is 



B 



f P 2 dV 
Jp%dV 



P{ P )^dp. 



(9) 



If the survival fraction /(pi) has a power-law depen- 
dence f(pi) oc Pi° , then there will be a contribution 
~ /(P©)(Pmax/po) 1 ~ a /(l — a) to the enhancement from 
the high-density tail in P(p). Thus, a large annihila- 
tion enhancement requires a < 1, /(p©) not too small 
(/ (P©) > 0.2), and p m ax 3> p. For example, we estimate 
that for f(p) = 0.2 (a constant) and p max ~ 1000 p©, 
the enhancement factor will be B ~ 200. The enhance- 
ment factor will be roughly proportional to /(p©), and 
it will decrease, possibly sharply, as a is increased from 
zero. Our numerical results find that B decreases to val- 
ues B < 10 for a = 0.5 and /(p ) < 0.2. It is also 
important to note that if the power-law exponent of a is 
such as to allow a large enhancement B, then the value 
of that enhancement is likely to depend strongly on the 
cutoff density p max . 

For the case where the Milky Way halo is composed 
of discrete subhalos with no smooth component, the an- 
nihilation rate in this halo is enhanced over that in the 
uniform-density halo by a factor 



B(c v ) = 



47r J r 2 dr [p(r)] 5 
47r J r 2 dr p 2 

1 Cyff(Cy) 

3[/(cv)] 2 ' 



(10) 



for c v ^> 1, in broad 



where g(c v ) = (1/3) [1 — (1 + c v )~ 3 ]. This enhancement 
factor is B(c v = 1) ~ 2.6, it grows to B(c v — 10) ~ 50, 
and grows roughly as c^/9[lnc v ] 2 
agreement with the results presented in Ref. [2l| 
In summary, we conclude that: 

(1) Very large (i.e., ^ 10) enhancements to the anni- 
hilation rate require a survival fraction / that is roughly 
constant with density. In other words, the survival frac- 
tion for a halo must be largely independent of its forma- 
tion time. If earlier halos are less likely to survive, then a 
large annihilation enhancement requires a survival frac- 
tion today that is large, and perhaps too large to be 
consistent with numerical simulations of halo formation. 
This is consistent with Ref. |2(| , who claim that survival 
fractions of the earliest halos are ~ 0.1% — 0.5% and an 
annihilation fraction B ~ 2 — 5. 

(2) The annihilation enhancement generally increases 
at the expense of the local density. For example, in the 
a = model, an annihilation enhancement B « 150 
implies a reduction by a factor of three in the local dark- 
matter density. The annihilation enhancement will gen- 
erally increase and the local smooth component will also 
generally decrease for larger p max . 



(3) Strictly speaking, the PDF derived here is for the 
dark-matter density locally (or at some other specified 
point), and the annihilation enhancement B is the en- 
hancement in the annihilation rate in a local volume. It 
is expected that the PDF, as well as B, will vary with 
radius in the Galactic halo. Since the central higher- 
density regions of the halo presumably formed earlier, 
the amount ot substructure should be reduced there (i.e., 
a larger smooth fraction and smaller B), and conversely 
for larger radii (smaller smooth fraction and larger B). 
These trends are consistent with the enhancement fac- 
tors found, for example, in Ref. [ijjj], which calculate the 
annihilation intensities as a function of observation di- 
rection. 



III. SIMULATION-INSPIRED RESULTS 

We now proceed to build on the analytic calculation 
in the previous Section by implementing a substructure 
fraction, a subhalo mass function, and a range of concen- 
tration parameters taken from simulations. 

We assume that when averaged over time, the solar 
neighborhood has a mean density p© ~ 0.4 GeVcm -3 . 
In a halo that contains substructure, this density can be 
split into a smooth component, and a component that 
arises from the presence of subhalos, 



P© 



Psm 
Psm 



Psub 
rn 



M - 



dn 



d\nM 



-dhiM. 



(11) 



Here, dn/d In M is the number density of objects per log- 
arithmic mass interval, and the integral is performed over 
a mass range, e < M < r\ < Mmw- Strictly speaking, 
the mass function depends on the primordial power spec- 
trum and upon the physics of halo merging and stripping 
of halos. 

Numerical simulations find that the mass function can 
be approximated over a wide mass rang e by a power law, 
dn/dhiM ~ M~@, with (3 ~ 0.9 [UTIM Hf, with in- 
dications that a similar form remains down to sub-solar 
subhalos [2d. |31( . We normalize the mass function so 
that a fraction £ of the local mean dark-matter density 
is in objects with masses between e w 10~ 10 M Q and 
i] w 10 12 M Q . This implies that the value of the smooth 
component is p sm = (1 — £)Po- The upper limit rj of 
the integration must always be less than the mass of the 
Milky Way. However, the lower integration limit e de- 
pends on the physics of the dark-matter particle and on 
small-scale structure assemblage. 

We model each subhalo with a two-parameter NFW 
profile, taking the two parameters to be the mass and 
the concentration parameter, and using a virial radius 
(maximum radius) defined by assuming the density of 
each halo is 200 times the matter density at the redshift 
of formation. Each subhalo is then assigned a concen- 
tration that is selected from a P(c v ) log-normal distribu- 
tion about a mean determined by the value of o~(M) and 
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the evolution of linear perturbations, and with a scat- 
ter <r[lpg(c v )] = 0.14 inferred from numerical simulations 
[32L [33l | . In this model, the mean concentration is a weak 
function of mass, c v ~ M~ 7 , where 7 w 0.13 for scales 
near M* at z sa 0, and scales as c v w 33(M/1O 8 M )~ - 06 
for halos with masses M < 10 8 M©. 

We calculate the density probability distribution func- 
tion as 



100 



PM-f 



dn dv(M,p) 



d\nM 



dp 



dlnM. 



(12) 



Here, dv(M, p)/dp = f[dv(M(c v ))/dp]P(c v )dc v , where 
v(M,p) — 47r(r s f c ) 3 /3 is the volume in a halo of mass 
M, where f c is obtained by solving 



r c (l + f c ) 2 = p s /p. 



(13) 



In solving Eq. (|13[) . we assume that halos can only have 
a finite size, given by f max — c v . If the solution is such 
that f < cv, then the volume with density greater than 
p is simply Airr^r 3 /3, while if f > c v , then the volume is 
given by 47rr 3 c 3 /3. However, in order to approximate the 
effects of tidal interactions, Eqs. (fl"2|) and (fl~3|) are always 
solved for p > p sm ; i.e., we assume the tidal radius of each 
subhalo is defined as the radius where the density of the 
halo is equal to the local smooth density component. 

Fig. [3] shows the density probability distribution func- 
tion as a function of p, expressed as a fraction of p@. 
Different line types depict the PDFs that corresponds to 
£ =0.1, 0.5, and 0.8. As before, we smooth the Dirac 
delta function at the value of the smooth density by a 
Gaussian with an RMS one-tenth the smooth compo- 
nent density. The range that corresponds to the sub- 
halo population for each £ represents the spread in P(p) 
that arises from a spread in the lower integration limit 
e and in the power-law exponent (3 of the subhalo mass 
function. More precisely, for each value of £, the upper 
curve corresponds to a subhalo power-law index (3 = 1.1 
and integration limits e = 10 -10 Mq and rj = 10~ 9 Mq; 
i.e. a steep subhalo mass function with most of the pop- 
ulation in objects of extremely small mass. The lower 
curve is obtained by flattening the subhalo mass func- 
tion to a power-law exponent f3 = 0.7 and integration 
limits e = 10 10 M Q (since we know the LMC exists) and 
7] = 10 12 M Q , roughly speaking the Milky Way mass; 
i.e., this is a model with only very massive substructure. 

The value of £ fixes the value of p sm , and densities 
p > Psm come from subhalos. The slope of P(p) is set in 
Eq. (fT2"]) by the slope of v(p). At densities where r <C 1, 
the volume within which the density is greater than a 



particular value is v 



For lower densities, 



probing the outer regions of a halo — i.e., f 3> 1 — the 



volume scales as v 



P 



Thus, we expect a slope 



for P(p) of —2 for low densities, and a slope of —4 as the 
density increases. 

The sharp transition between the smooth component 
and the subhalos in Fig. [3] results from the assumption 
of a pure power law for the dark-matter profile in sub- 
halos Q. In reality, the transition will be much more 



3 




0.001 f 



0.0001 



100 



FIG. 3: The probability distribution function in the solar 
neighborhood for the simulation-inspired calculation. The 
solid, long-dash and dot-dash curves correspond to £ = 0.1, 
0.5, and 0.8 respectively. As in Fig. [TJ we smooth the Dirac 
delta function for the density value of the smooth component 
with a Gaussian of rms a tenth of the smooth-component den- 
sity. The upper and lower curves for each value of £ show the 
range of contribution of the subhalo population that arises 
from uncertainties in the subhalo mass function, as well as 
the subhalo population (see text). 



gradual, reflecting the fact that subhalos are embedded 
in the potential well of their host Milky Way halo. As 
£ — » 1, the smooth component approaches zero. There 
is still a lower limit, p > 0.04 GeVcm~ 3 , attained at a 
value £ > 0.92, to the local halo density that arises from 
the overlap of the outskirts of subhalos. These calcula- 
tions thus imply that there is significant probability for 
the local halo density to be below p Q , but also that there 
is likely a minimum possible value to the local density. 
Simulations (e.g., Ref. [ll[) suggest a value £ ~ 0.5, but 
further investigation is required to determine this impor- 
tant parameter more precisely and robustly. 

The annihilation enhancement in this approach can be 
estimate as in Ref. [2l[ . Namely, in a halo of mass M and 
with a cut-off in the subhalo mass function at mo, the 
boost factor is the solution of an integral equation that 
takes into account the halo-in-halo problem, and is ap- 
proximated by B « 0.1[(M/mo) 013 — 1]. For a cut-off in 
the power spectrum at microhalo scales, mo « 1O _6 M0, 
and a Milky Way mass of M ps 1O 12 M , B ps 20. The 
weak dependence of the boost factor to the cut-off scale 
of the subhalo mass function is an outcome of the flatness 
of the dark-matter power spectrum, which manifests it- 
self in the concentration-mass relationship which enters 
the boost-factor calculation (see Ref. (2lj). 
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IV. DISCUSSION 

There is a vast literature on galactic substructure, 
much of it discussing the implications for dark-matter 
detection. Most of this work focuses on the annihilation 
boost factor, but there are few papers that discuss the 
implications for direct detection. The primary focus of 
these papers is then on the effects of substructure on the 
local dark-matter velocity distribution, with less atten- 
tion to the possible implications for the local dark-matter 
density. Here we review some of the work that discusses 
the implications of substructure on the local dark-matter 
density. 

In Ref. [231 ]. the authors use numerical simulations to 
investigate the effects of substructure on the local veloc- 
ity distribution and on the local density. The conclusion 
of that paper is that most of the local dark matter is 
smoothly distributed (illustrated in their Fig. 5, which 
resolves structure down to ~ 1O 7 M0). The fraction of 
mass that is bound in halos less massive than 10 7 M Q is 
less than 4%; this is the fraction of the mass that could 
have survived the tidal field of the Milky Way. They thus 
conclude that it is highly unlikely that a good fraction of 
the halo is in Earth-mass subhalos. However, their calcu- 
lation does not consider the halos-in-halos problem: i.c, 
their mass function dN / dm keeps track only of the most 
massive halo in which a given particle resides, not the 
least massive subhalo. Put another way, their calculation 
counts only the Earth-mass halos that did not get incor- 
porated into more massive halos, and it disregards those 
that did, assuming simply that they were completely dis- 
rupted. Still, this paper does argue that substructures 
at later stages in the merger hierarchy get largely erased. 
Although isolated Earth-mass substructures may survive 
in a smooth Milky Way halo (as pointed out also in 
Ref. (22| ) , when scaled to earlier generations in the hierar- 
chy, the conclusions of this paper imply that Earth-mass 
subhalos are likely to be disrupted when they merge into 
subsequent stages in the hierarchy (e.g., into 1000 Mq 
halos). If so, then f(p) will be very small (f(p) <C 0.2)). 

A local dark matter density PDF P(p) is calculated 
in Ref. [25j], but it is a different distribution. In partic- 
ular, their calculation assumes that the vast majority of 
the local dark matter is smoothly distributed, and that 
only 1% — 5% of the local dark matter may be in sub- 
structure. The density p that appears in their PDF is 
thus the density of this additional substructure compo- 
nent, which they assume comes from the latest subhalos 
accreted onto the Milky Way halo. Their mass function 
does extend to small masses, but they do not consider 
halos-in-halos. Ref. [26| is similar in spirit, but consid- 
ers in particular the effects of the tails of the Sagittarius 
dwarf. 

The approach that most closely resembles ours is per- 
haps that in Ref. [22j • They recognize the scale- invariant 
nature of the problem — that is, that there may be ha- 
los in halos — but then note that the resolution limits of 
their simulation prohibits them from making definitive 



claims about earth-mass objects. Their calculations in- 
dicate that the singular cores of subhalos may always 
survive, even if most of the mass from a given subhalo is 
stripped. They provide as an example a subhalo orbiting 
at 20 kpc in the Milky Way in which only 0.3% of the 
initial mass remains after four orbits, but then show that 
the survival fraction may be as high as 40% for an orbit at 
40 kpc. The simulations of Ref. [22j make the important 
point that tidal tails broaden rapidly, and this justifies 
the assumption of our first analytic model that matter 
stripped from a subhalo is rapidly smoothly distributed 
in the new larger halo. 

This assumption is also justified by the recent results 
presented in Ref. (24[. In this work, the authors pro- 
posed a technique for calculating the fine-grained phase 
space structure in dark matter halos from cosmological 
N-body simulations. In demonstrating the effectiveness 
of this new method, they studied the evolution of the 
dark matter density that arises from streams in NFW- 
like potentials, and found that a very large number of 
streams (~ 10 5 ) may potentially be present in the solar 
neighborhood. If each stream has a diameter of order 
^kpc, it means that their entanglement and evolution 
has an effect in the local smooth dark matter density, 
and thus do not address the potential enhancement on 
much smaller scales. 



V. CONCLUSIONS AND CHALLENGES 

It is evident from simulations and analytic arguments 
that some fraction of the local Milky Way dark mat- 
ter may be in subhalos. The implications of substruc- 
ture for indirect detection of WIMPs have been studied 
broadly, with the conclusion that there may be large en- 
hancements in annihilation rates over the rates predicted 
assuming a smooth halo. The implications for direct 
searches have, however, been largely overlooked. This 
is a possibly serious omission, as one consequence of sub- 
structure is that the local density will be smaller than 
the smoothly-distributed local density usually assumed. 

We have taken a few first steps to understand the im- 
plications of substructure for the local density. Our cen- 
tral goal is a calculation of the PDF P{p) for the local 
dark-matter density p. This P(p) will need to be taken 
into account when interpreting the implications of null 
dark-matter searches for constraints to the particle-dark- 
matter parameter space (e.g., couplings and/or elastic- 
scattering cross sections). 

We considered two simple scenarios for substructure: 
In the first, early generations of very dense subhalos sur- 
vive with some probability. The advantage of this ap- 
proach is that if subhalo survival fractions can be mea- 
sured in simulations for recently merged subhalos, the 
results might be extrapolated to the much earlier gener- 
ations (much smaller and denser subhalos) that may be 
below the resolution of simulations. 

In the second approach, the halo is assumed to consist 
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of recently formed subhalos, each with an NFW profile. 
This approach provides a simple, albeit approximate, 
way to understand the effects of larger-scale substruc- 
ture, from more recent stages in the merger hierarchy, 
on the PDF. We then pursued this approach further us- 
ing subhalo mass functions and concentration-parameter 
distributions taken from simulations. This approach does 
not take into account the possible contribution of subha- 
los within subhalos. In principle, a complete solution 
for the PDF, including substructures on the largest and 
smallest scales, can be obtained by convolving our two 
calculations. This, however, will be left for future re- 
search. 

Substructure scenarios that yield larger annihilation 
enhancements generally imply a smaller local dark- 
matter density. Very large annihilation enhancements 
require that the very densest substructures, which gen- 
erally form earlier, must survive through all later gen- 
erations of structure formation. If earlier substructures 
are less likely to survive than more recent substructures, 
then a very large annihilation enhancement is unlikely. 

So, how small can the local density be? The smallest 
local density in the models we surveyed was one tenth the 
canonical value of 0.4 GeVcm -3 , usually assumed for a 
smoothly distributed halo. This small value was obtained 
from our simulation-inspired result using what we be- 
lieve to be an overly conservative estimate of the smooth 
fraction, which we obtained by truncating NFW subha- 
los when their density falls below the mean halo density. 
More realistic values for the smooth fraction are probably 
in the range of 50%-80%, which would then correspond in 
our simulation-inspired results to the £ = 0.8 and £ = 0.5 
distributions, respectively, shown in Fig. [31 The smallest 
local density in our analytic model for early substructure 
was 0.3, obtained using the high value, /(p©) = 0.2, and 
assuming a constant /(pi). If, however, /(p©) is lower, 
then the local density is increased. More importantly, 
it is quite likely that f(p) decreases with p, and if so, 
the local density is increased to ~ 0.8 times its canonical 
value, even for /(p©) = 0.2, and even higher for smaller 
/(p©)- A combination of our two approaches, to take 
into account both low- and high-density substructures, 
will likely show that the local density is no less than half 
the canonical value. 

The volume of the halo probed during a three-year 
direct-detection experiment is very small, and so the halo 
density p that we have been discussing can be safely as- 
sumed to be the density averaged over the duration of a 
direct-detection experiment. However, the rate for pro- 
duction of energetic neutrinos from WIMP annihilation 
in the Sun or Earth lags behind the rate for capture of 
WIMPs from the halo by an equilibration time t eq [34l[35| 
that can vary considerably with the WIMP's mass and 
clastic-scattering and annihilation cross sections; typical 
values might be t® ~ 5 X 10 7 yr and t® - 10 10 yr 
An energetic-neutrino search thus probes the halo den- 
sity averaged over a much larger volume (for the latest 
bounds on the flux of high-energy neutrinos from anni- 



hilations in the Earth and the Sun see [33, [H, HH). The 
halo density averaged over this volume will thus have a 
PDF that will be narrower than that for direct detection, 
and the minimum density will be closer to the mean halo 
density. However, even though initially it was assumed 
that the velocity distribution function of dark matter par- 
ticles mirrors that in free space (thus justifying the use 
of a Maxwell-Boltzmann distribution) [40j], it is possible 
that the velocity distribution function is surpressed in the 
low-velocity tail due to the effects of solar capture and 
WIMP diffusion in the solar system due to the presence 
of other planets [4l], |42j • This surpression manifests it- 
self as a reduction in the number density of dark matter 
particles near the Earth for WIMPs with masses greater 
than few hundred GeV. Nevertheless, there is a possi- 
bility that for low-mass WIMPs, given the longer equili- 
bration time for the Earth relative to that for the Sun, 
the annihilation signal from the Earth could be boosted 
relative to that for the Sun, if the Solar System passed 
through a very dense subhalo at a time t® < t < tf 
ago. We leave a more detailed calculation of the density 
smoothed over timescales relevant for energetic-neutrino 
searches to future work (43|. 

The PDF P(p) allows us to evaluate also the boost 
factor in the annihilation rate. It must be emphasized, 
though, that the boost factor B we have considered is 
the boost only in the local annihilation rate (per unit 
volume). Strictly, speaking, the PDF may vary from one 
position in the halo to another. Thus, for example, if we 
assume that the halo is spherically symmetric (after aver- 
aging over substructure fluctuations), then the PDF will 
be a function P(p; r) of Galactocentric radius r, as well 
as the density p. The mean density of this distribution, 
as a function of r, will be the density of the smooth NFW 
profile that best fits the rotation curve. Qualitatively, we 
expect that the high-density tail will be more pronounced 
(more substructure) at larger r and less pronounced (less 
substructure) at smaller r. The PDF P(p; r) can then 
be used to determine the boost factor as a function of 
position in the halo; again, we expect B(r) to vary with 
r and to generally increase with r. This B(r) will be es- 
sential to compute annihilation fluxes along various lines 
of sight through the Galaxy. We encourage future au- 
thors to describe the effects of substructure in terms of 
P{p\ r), or at least in terms of B(r), as this may facilitate 
comparison between the conclusions of different studies. 
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